Post-COMPASS Plots
Malisa Smith
version June 16, 2020
Consider the cohort as a single group to start. Do not stratify by case-control status.
Stack the heatmaps for CD4. Are there obvious patterns to the antigens?
Stack the heatmaps for NOTCD4. Are there patterns? Are these patterns different from CD4?
Calculate FS/PFS for each stimulation in CD4 and compare side-by-side in boxplots and test by ANOVA. Is there a clear winner? Calculate FS/PFS for each stimulation in NOTCD4 and compare with boxplot and ANOVA. Is there a winner? Is this different than CD4?
It may also make sense to make a corrplot of the FS/PFS from the different stims. This is asking a different question though, i.e. whether a patient’s cytokine responses tend to be consistently elevated/not-elevated across multiple stims.
Regress FS/PFS against age (linear). Is there an effect of age on T cell functionality? Regress FS/PFS against days from symptom onset (linear). Is there a decrease in functionality with time? Stratify FS/PFS by sex and test by t-test (or Mann-Whitney, though sample size is large enough). Is there a difference?
Only then would you start asking the case-control questions.
library(here)
library(tidyverse)
library(COMPASS)
library(ggplot2)
library(data.table)
library(broom)
library(psych)
library(corrplot)
library(ggpubr)
library(knitr)
library(kableExtra)
library(cowplot)
source(here::here("scripts/my_pheatmap_and_plotMeanGamma.R"))
source(here::here("scripts/20200604_Helper_Functions.R"))save_output <- FALSERead in data
# Next 4 lines from Run_COMPASS.R
stims_for_compass_runs <- c("VEMP", "Spike 1", "Spike 2", "NCAP")
parent_nodes_for_compass_runs <- c("4+", "NOT4+")
stims_for_compass_runs_rep <- rep(stims_for_compass_runs, each = length(parent_nodes_for_compass_runs))
parent_nodes_for_compass_runs_rep <- rep(parent_nodes_for_compass_runs, times = length(stims_for_compass_runs))
crList <- purrr::pmap(.l = list(parent_nodes_for_compass_runs_rep,
stims_for_compass_runs_rep),
.f = function(parent, currentStim) {
currentStimForFile <- gsub(" ", "_", currentStim)
crPath <- here::here(sprintf("out/CompassOutput/%s/%s/COMPASSResult_%s_%s.rds", parent, currentStimForFile, parent, currentStimForFile))
readRDS(crPath)
}) %>%
set_names(gsub("+", "", gsub(" ", "_", paste0(parent_nodes_for_compass_runs_rep, "_", stims_for_compass_runs_rep)), fixed=TRUE))
names(crList)## [1] "4_VEMP" "NOT4_VEMP" "4_Spike_1" "NOT4_Spike_1" "4_Spike_2"
## [6] "NOT4_Spike_2" "4_NCAP" "NOT4_NCAP"CD4RunNames <- c("4_Spike_1", "4_Spike_2", "4_NCAP", "4_VEMP")
NotCD4RunNames <- c("NOT4_Spike_1", "NOT4_Spike_2", "NOT4_NCAP", "NOT4_VEMP")Stacked heatmaps
plotCOMPASSResultStack(crList[CD4RunNames],
row_annotation = c("Stim"),
variable = "Stim",
show_rownames = FALSE,
main = "CD4 COMPASS Runs",
fontsize = 15, fontsize_row = 15, fontsize_col = 13)## The 'threshold' filter has removed 6 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&IL17a Ax700&!IL2 PE&!IL4/5/13 APC&!TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&!TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&IL4/5/13 APC&!TNFa FITC
plotCOMPASSResultStack(crList[NotCD4RunNames],
row_annotation = c("Stim"),
variable = "Stim",
show_rownames = FALSE,
main = "Not-CD4 COMPASS Runs",
fontsize = 15, fontsize_row = 15, fontsize_col = 13)## The 'threshold' filter has removed 2 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&IL17a Ax700&!IL2 PE&IL4/5/13 APC&!TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC
plotCOMPASSResultStack(crList[c(CD4RunNames, NotCD4RunNames)],
row_annotation = c("Stim"),
variable = "Stim",
show_rownames = FALSE,
main = "CD4 and Not-CD4 COMPASS Runs",
fontsize = 15, fontsize_row = 15, fontsize_col = 13)## The 'threshold' filter has removed 3 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&IL4/5/13 APC&!TNFa FITC
if(save_output) {
png(filename=here::here("out/PostCompassPlots/20200613_CD4_COMPASS_Heatmap_All_Stims.png"),
width=700, height=1000)
plotCOMPASSResultStack(crList[CD4RunNames],
row_annotation = c("Stim"),
variable = "Stim",
show_rownames = FALSE,
main = "CD4 COMPASS Runs",
fontsize = 14, fontsize_row = 13, fontsize_col = 13)
dev.off()
png(filename=here::here("out/PostCompassPlots/20200613_NotCD4_COMPASS_Heatmap_All_Stims.png"),
width=700, height=1000)
plotCOMPASSResultStack(crList[NotCD4RunNames],
row_annotation = c("Stim"),
variable = "Stim",
show_rownames = FALSE,
main = "Not-CD4 COMPASS Runs",
fontsize = 14, fontsize_row = 13, fontsize_col = 13)
dev.off()
png(filename=here::here("out/PostCompassPlots/20200613_CD4_and_NotCD4_COMPASS_Heatmap_All_Stims.png"),
width=900, height=2000)
plotCOMPASSResultStack(crList[c(CD4RunNames, NotCD4RunNames)],
row_annotation = c("Stim"),
variable = "Stim",
show_rownames = FALSE,
main = "CD4 and Not-CD4 COMPASS Runs",
fontsize = 14, fontsize_row = 13, fontsize_col = 13)
dev.off()
}FS/PFS across Stims
fs_pfs_df <- lapply(c(CD4RunNames, NotCD4RunNames), function(n) {
as.data.frame(COMPASS::FunctionalityScore(crList[[n]])) %>%
rownames_to_column("SAMPLE ID") %>%
rename(FS = `COMPASS::FunctionalityScore(crList[[n]])`) %>%
left_join(as.data.frame(COMPASS::PolyfunctionalityScore(crList[[n]])) %>%
rownames_to_column("SAMPLE ID") %>%
rename(PFS = `COMPASS::PolyfunctionalityScore(crList[[n]])`),
by = "SAMPLE ID") %>%
mutate(COMPASS_run = n)
} ) %>%
data.table::rbindlist() %>%
separate(COMPASS_run, into = c("parent", "STIM"), remove = F, extra = "merge") %>%
mutate(STIM = factor(STIM, levels = c("Spike_1", "Spike_2", "NCAP", "VEMP")))#,
#`SAMPLE ID` = factor(`SAMPLE ID`, levels = unique(`SAMPLE ID`)))
table(fs_pfs_df$STIM, fs_pfs_df$`SAMPLE ID`) # remove 07B for quade.test in order to have complete block design##
## 07B 102C 103C 106C 109C 112C 113C 117C 121C 122C 127C 128C 12C 13
## Spike_1 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## Spike_2 0 2 2 2 2 2 2 2 2 2 2 2 2 2
## NCAP 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## VEMP 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##
## 130C 131C 133C 136C 139C 141C 142C 143C 147C 148C 150C 15514 15518
## Spike_1 2 2 2 2 2 2 2 2 2 2 2 2 2
## Spike_2 2 2 2 2 2 2 2 2 2 2 2 2 2
## NCAP 2 2 2 2 2 2 2 2 2 2 2 2 2
## VEMP 2 2 2 2 2 2 2 2 2 2 2 2 2
##
## 15529 15530 15545 15548 15554 15575 157C 159C 161C 162C 164C 169C
## Spike_1 2 2 2 2 2 2 2 2 2 2 2 2
## Spike_2 2 2 2 2 2 2 2 2 2 2 2 2
## NCAP 2 2 2 2 2 2 2 2 2 2 2 2
## VEMP 2 2 2 2 2 2 2 2 2 2 2 2
##
## 170C 1C 21C 23 25 25C 29C 36C 38C 51C 56C 59C 69C 6C 76C 80C 90C 93C
## Spike_1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## Spike_2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## NCAP 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## VEMP 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##
## 96C BWT20 CLT1 HS09 HS10 HS8
## Spike_1 2 2 2 2 2 2
## Spike_2 2 2 2 2 2 2
## NCAP 2 2 2 2 2 2
## VEMP 2 2 2 2 2 2fs_pfs_df.complete <- fs_pfs_df %>%
dplyr::filter(`SAMPLE ID` != "07B") %>%
mutate(`SAMPLE ID` = factor(`SAMPLE ID`, levels = unique(`SAMPLE ID`))) %>%
arrange(STIM, `SAMPLE ID`)
fs_pfs_df.complete.4 <- fs_pfs_df.complete %>% dplyr::filter(parent == "4")
fs_pfs_df.complete.not4 <- fs_pfs_df.complete %>% dplyr::filter(parent == "NOT4")
# Perform statistical tests for Functionality Score across Stims
# Quade test is to wilcoxon signed rank test like ANOVA is to t-test
quade_fs_4 <- quade.test(FS ~ STIM | `SAMPLE ID`, data = fs_pfs_df.complete.4)
quade_fs_4##
## Quade test
##
## data: FS and STIM and SAMPLE ID
## Quade F = 73.3, num df = 3, denom df = 183, p-value < 2.2e-16quade_fs_not4 <- quade.test(FS ~ STIM | `SAMPLE ID`, data = fs_pfs_df.complete.not4)
quade_fs_not4##
## Quade test
##
## data: FS and STIM and SAMPLE ID
## Quade F = 41.83, num df = 3, denom df = 183, p-value < 2.2e-16pw_fs_4 <- pairwise.wilcox.test(fs_pfs_df.complete.4$FS,
fs_pfs_df.complete.4$STIM,
p.adjust.method="bonferroni",
paired=TRUE)
pw_fs_4$p.value## Spike_1 Spike_2 NCAP
## Spike_2 1.632616e-01 NA NA
## NCAP 4.082570e-02 3.720152e-06 NA
## VEMP 6.244948e-11 4.659364e-11 6.244948e-11broom::tidy(pw_fs_4)## # A tibble: 6 x 3
## group1 group2 p.value
## <chr> <chr> <dbl>
## 1 Spike_2 Spike_1 1.63e- 1
## 2 NCAP Spike_1 4.08e- 2
## 3 VEMP Spike_1 6.24e-11
## 4 NCAP Spike_2 3.72e- 6
## 5 VEMP Spike_2 4.66e-11
## 6 VEMP NCAP 6.24e-11pw_fs_not4 <- pairwise.wilcox.test(fs_pfs_df.complete.not4$FS,
fs_pfs_df.complete.not4$STIM,
p.adjust.method="bonferroni",
paired=TRUE)
pw_fs_not4$p.value## Spike_1 Spike_2 NCAP
## Spike_2 1.000000e+00 NA NA
## NCAP 7.599030e-09 1.461474e-09 NA
## VEMP 1.984327e-02 2.881666e-01 4.411433e-10broom::tidy(pw_fs_not4)## # A tibble: 6 x 3
## group1 group2 p.value
## <chr> <chr> <dbl>
## 1 Spike_2 Spike_1 1.00e+ 0
## 2 NCAP Spike_1 7.60e- 9
## 3 VEMP Spike_1 1.98e- 2
## 4 NCAP Spike_2 1.46e- 9
## 5 VEMP Spike_2 2.88e- 1
## 6 VEMP NCAP 4.41e-10# Perform statistical tests for PolyFunctionality Score across Stims
quade_pfs_4 <- quade.test(PFS ~ STIM | `SAMPLE ID`, data = fs_pfs_df.complete.4)
quade_pfs_4##
## Quade test
##
## data: PFS and STIM and SAMPLE ID
## Quade F = 63.947, num df = 3, denom df = 183, p-value < 2.2e-16quade_pfs_not4 <- quade.test(PFS ~ STIM | `SAMPLE ID`, data = fs_pfs_df.complete.not4)
quade_pfs_not4##
## Quade test
##
## data: PFS and STIM and SAMPLE ID
## Quade F = 32.722, num df = 3, denom df = 183, p-value < 2.2e-16pw_pfs_4 <- pairwise.wilcox.test(fs_pfs_df.complete.4$PFS,
fs_pfs_df.complete.4$STIM,
p.adjust.method="bonferroni",
paired=TRUE)
pw_pfs_4$p.value## Spike_1 Spike_2 NCAP
## Spike_2 3.831211e-02 NA NA
## NCAP 1.000000e+00 6.717394e-04 NA
## VEMP 7.960767e-11 4.659364e-11 6.556256e-11# broom::tidy(pw_pfs_4)
pw_pfs_not4 <- pairwise.wilcox.test(fs_pfs_df.complete.not4$PFS,
fs_pfs_df.complete.not4$STIM,
p.adjust.method="bonferroni",
paired=TRUE)
pw_pfs_not4$p.value## Spike_1 Spike_2 NCAP
## Spike_2 1.000000e+00 NA NA
## NCAP 2.137504e-08 1.518882e-08 NA
## VEMP 1.000000e+00 1.000000e+00 3.142649e-09# broom::tidy(pw_pfs_not4)
# Now plot the results
# Note that the dataset used in statistics excludes 07B, but 07B is plotted below (and the median bar is calculated including 07B)
parent.labs <- c("CD4", "Not-CD4")
names(parent.labs) <- c("4", "NOT4")
plot_fs_pfs_vs_stim <- function(fs_or_pfs) {
ggplot(fs_pfs_df,
aes(STIM, !!as.name(fs_or_pfs))) +
theme_bw(base_size = 22) +
geom_line(aes(group = `SAMPLE ID`), color="grey") +
geom_point(color="grey") +
facet_grid(. ~ parent,
labeller = labeller(parent = parent.labs)) +
geom_errorbarh(data = fs_pfs_df %>% group_by(parent, STIM) %>% summarise(med = median(!!as.name(fs_or_pfs))),
aes(y = med,
xmax = 1.5 + 0.35,
xmin = 1.5 - 0.35,
height = 0),
position=position_dodge(width=0), color = "black") +
labs(title = sprintf("%s vs Stim", fs_or_pfs),
y = fs_or_pfs) +
theme(axis.title.x = element_blank(),
axis.title.y = element_text(size=22),
axis.text.y = element_text(color="black", size=22),
axis.text.x = element_text(color="black", size=22),
plot.title = element_text(hjust = 0.5),
panel.grid = element_blank(),
legend.position = "none",
plot.margin = margin(1.3, 0.2, 0, 0.2, "cm")) +
scale_x_discrete(labels=c(c("Spike_1" = "Spike 1", "Spike_2" = "Spike 2")),
expand = c(0.1,0.1))
}
plot_fs_pfs_vs_stim("FS")
# Remembering these p-values are adjusted
as.data.frame(pw_fs_4$p.value) %>%
rownames_to_column('tmp') %>%
mutate_all(~cell_spec(.x, color = ifelse(!is.na(.x) & .x < 0.05, "red", "black"))) %>%
column_to_rownames('tmp') %>%
kable(format = "html", escape = F, row.names = T) %>%
kable_styling("striped", full_width = F)| Spike_1 | Spike_2 | NCAP | |
|---|---|---|---|
| Spike_2 | 0.163261630177697 | NA | NA |
| NCAP | 0.0408257009089326 | 3.72015208085264e-06 | NA |
| VEMP | 6.24494803673655e-11 | 4.65936369334732e-11 | 6.24494803673655e-11 |
as.data.frame(pw_fs_not4$p.value) %>%
rownames_to_column('tmp') %>%
mutate_all(~cell_spec(.x, color = ifelse(!is.na(.x) & .x < 0.05, "red", "black"))) %>%
column_to_rownames('tmp') %>%
kable(format = "html", escape = F, row.names = T) %>%
kable_styling("striped", full_width = F)| Spike_1 | Spike_2 | NCAP | |
|---|---|---|---|
| Spike_2 | 1 | NA | NA |
| NCAP | 7.5990300522029e-09 | 1.46147350795162e-09 | NA |
| VEMP | 0.0198432725984609 | 0.288166634062905 | 4.41143306237813e-10 |
plot_fs_pfs_vs_stim("PFS")
as.data.frame(pw_pfs_4$p.value) %>%
rownames_to_column('tmp') %>%
mutate_all(~cell_spec(.x, color = ifelse(!is.na(.x) & .x < 0.05, "red", "black"))) %>%
column_to_rownames('tmp') %>%
kable(format = "html", escape = F, row.names = T) %>%
kable_styling("striped", full_width = F)| Spike_1 | Spike_2 | NCAP | |
|---|---|---|---|
| Spike_2 | 0.0383121094257177 | NA | NA |
| NCAP | 1 | 0.000671739390258177 | NA |
| VEMP | 7.96076669051054e-11 | 4.65936369334732e-11 | 6.55625616466633e-11 |
as.data.frame(pw_pfs_not4$p.value) %>%
rownames_to_column('tmp') %>%
mutate_all(~cell_spec(.x, color = ifelse(!is.na(.x) & .x < 0.05, "red", "black"))) %>%
column_to_rownames('tmp') %>%
kable(format = "html", escape = F, row.names = T) %>%
kable_styling("striped", full_width = F)| Spike_1 | Spike_2 | NCAP | |
|---|---|---|---|
| Spike_2 | 1 | NA | NA |
| NCAP | 2.13750379422692e-08 | 1.51888161479098e-08 | NA |
| VEMP | 1 | 1 | 3.14264895925017e-09 |
As we saw in the heatmaps, NAP and Spike 1 consistently have the highest FS/PFS scores, and VEMP has the lowest.
How correlated are FS/PFS across stims?
fs_pfs_df_wide <- fs_pfs_df %>%
pivot_wider(id_cols = c(`SAMPLE ID`),
names_from = c(STIM, parent),
values_from = c(FS, PFS))
pairs.panels(fs_pfs_df_wide %>%
dplyr::select(-c(`SAMPLE ID`)) %>%
dplyr::select(matches("^FS.*_4")),
method = "pearson", # correlation method, "r" displayed in upper off diagonals. Not r^2.
hist.col = "#00AFBB",
density = TRUE, # show density plots,
lm = T,
main = "CD4 FS Correlation across Stims")
pairs.panels(fs_pfs_df_wide %>%
dplyr::select(-c(`SAMPLE ID`)) %>%
dplyr::select(matches("^FS.*_NOT4")),
method = "pearson", # correlation method, "r" displayed in upper off diagonals. Not r^2.
hist.col = "#00AFBB",
density = TRUE, # show density plots,
lm = T,
main = "Not-CD4 FS Correlation across Stims")
pairs.panels(fs_pfs_df_wide %>%
dplyr::select(-c(`SAMPLE ID`)) %>%
dplyr::select(matches("^PFS.*_4")),
method = "pearson", # correlation method, "r" displayed in upper off diagonals. Not r^2.
hist.col = "#00AFBB",
density = TRUE, # show density plots,
lm = T,
main = "CD4 FS Correlation across Stims")
pairs.panels(fs_pfs_df_wide %>%
dplyr::select(-c(`SAMPLE ID`)) %>%
dplyr::select(matches("^PFS.*_NOT4")),
method = "pearson", # correlation method, "r" displayed in upper off diagonals. Not r^2.
hist.col = "#00AFBB",
density = TRUE, # show density plots,
lm = T,
main = "Not-CD4 PFS Correlation across Stims")
M <- cor(fs_pfs_df_wide %>%
dplyr::select(-c(`SAMPLE ID`)) %>%
dplyr::select(matches("PFS")),
use = "complete.obs") # missing values are handled by casewise deletion
# mat : is a matrix of data
# ... : further arguments to pass to the native R cor.test function
cor.mtest <- function(mat, ...) {
mat <- as.matrix(mat)
n <- ncol(mat)
p.mat<- matrix(NA, n, n)
diag(p.mat) <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
tmp <- cor.test(mat[, i], mat[, j], ...)
p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
}
}
colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
p.mat
}
getOption("na.action") # Default behavior of cor.test## [1] "na.omit"# matrix of the p-value of the correlation
p.mat <- cor.mtest(fs_pfs_df_wide %>%
dplyr::select(-c(`SAMPLE ID`)) %>%
dplyr::select(matches("PFS")))
head(p.mat[, 1:5])## PFS_Spike_1_4 PFS_Spike_2_4 PFS_NCAP_4 PFS_VEMP_4
## PFS_Spike_1_4 0.000000e+00 3.557519e-07 1.252908e-08 0.0003300951
## PFS_Spike_2_4 3.557519e-07 0.000000e+00 1.010421e-05 0.0196248181
## PFS_NCAP_4 1.252908e-08 1.010421e-05 0.000000e+00 0.0018545939
## PFS_VEMP_4 3.300951e-04 1.962482e-02 1.854594e-03 0.0000000000
## PFS_Spike_1_NOT4 2.113868e-09 4.275608e-03 6.103202e-03 0.0173764534
## PFS_Spike_2_NOT4 2.010616e-01 5.250470e-04 2.707909e-01 0.1543216526
## PFS_Spike_1_NOT4
## PFS_Spike_1_4 2.113868e-09
## PFS_Spike_2_4 4.275608e-03
## PFS_NCAP_4 6.103202e-03
## PFS_VEMP_4 1.737645e-02
## PFS_Spike_1_NOT4 0.000000e+00
## PFS_Spike_2_NOT4 2.657072e-01col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(M,
method="color", col=col(200),
type="upper", order="hclust",
addCoef.col = "black", # Add coefficient of correlation
tl.col="black", tl.srt=45, #Text label color and rotation
# Combine with significance
p.mat = p.mat, sig.level = 0.01, insig = "blank",
# hide correlation coefficient on the principal diagonal
diag=FALSE)
# significant coefficients are coloredObservations:
(1) Note that FS and PFS scores for the same condition tend to be highly correlated (not shown). Not surprising.
(2) Within PFS only, there is a mixed bag of significant associations occuring “between CD4 and Not-CD4 for the same condition” and “between different stims”. Tends to be much lower r than those for (1) though.
In general, we can’t rely on a particular condition to predict the FS/PFS of all the other conditions. Each STIM and co-receptor subset provides different information. (If the former were true, it might have simplified our upcoming analyses)
Regress FS/PFS vs Demographics
# Read in the final patient manifest with more complete data
data_manifest <- read.csv(here::here("data/Seshadri_HAARVI_PBMC_manifest_merged_11June2020.csv"), check.names = F, stringsAsFactors = F)
setdiff(sort(unique(fs_pfs_df$`SAMPLE ID`)), sort(unique(data_manifest$`Record ID`)))## [1] "07B" "BWT20" "CLT1" "HS09"setdiff(sort(unique(data_manifest$`Record ID`)), sort(unique(fs_pfs_df$`SAMPLE ID`)))## [1] "116C" "37C" "551432" "75C" "HS9"# Remember 116C and 37C didn't get run through COMPASS
# 75C was not run at all
# 07B/551432 didn't get run for Spike 2
fs_pfs_df_with_manifest <- fs_pfs_df %>%
mutate("Record ID" = dplyr::recode(`SAMPLE ID`,
"07B" = "551432",
"HS09" = "HS9")) %>% # also, 75C did not get run
full_join(data_manifest, by = c("Record ID" = "Record ID")) %>%
dplyr::filter(!(`Record ID` %in% c("116C", "37C", "75C"))) %>%
mutate(Days_Symptom_Onset_to_Visit_1 = as.numeric(`Days symptom onset to visit 1`),
Cohort = factor(ifelse(Cohort %in%
c(NA, "Healthy control", "Healthy control 2017-2018"), "Healthy", Cohort),
levels = c("Healthy", "Non-hospitalized", "Hospitalized")))## Warning: NAs introduced by coercionfs_pfs_df_with_manifest %>% dplyr::filter(is.na(FS) | is.na(`Record ID`) | is.na(`SAMPLE ID`) | is.na(STIM) | is.na(Cohort))## [1] SAMPLE ID FS
## [3] PFS COMPASS_run
## [5] parent STIM
## [7] Record ID Sample ID
## [9] Collection date Cell count
## [11] Cohort Age
## [13] Sex Race
## [15] Hispanic? Days symptom onset to visit 1
## [17] Pair ID Race_v2
## [19] Days_Symptom_Onset_to_Visit_1
## <0 rows> (or 0-length row.names)FS/PFS vs Age
# Note that p-values on ggscatter plots below are *not* adjusted, e.g.:
broom::tidy(lm(FS ~ Age, data = fs_pfs_df_with_manifest %>% dplyr::filter(parent == "4" & STIM == "Spike_2")))## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.0397 0.00704 5.64 0.000000529
## 2 Age 0.000230 0.000128 1.80 0.0772ggscatter(fs_pfs_df_with_manifest,
x = "Age", y = "FS",
color = "STIM", palette = "jco",
add = "reg.line") +
facet_grid(parent ~ STIM) +
stat_cor()## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 16 rows containing non-finite values (stat_smooth).## Warning: Removed 16 rows containing non-finite values (stat_cor).## Warning: Removed 16 rows containing missing values (geom_point).
ggscatter(fs_pfs_df_with_manifest,
x = "Age", y = "PFS",
color = "STIM", palette = "jco",
add = "reg.line") +
facet_grid(parent ~ STIM) +
stat_cor()## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 16 rows containing non-finite values (stat_smooth).## Warning: Removed 16 rows containing non-finite values (stat_cor).## Warning: Removed 16 rows containing missing values (geom_point).
FS/PFS vs Days from Symptom Onset to Visit
# Note that p-values on ggscatter plots below are *not* adjusted, e.g.:
broom::tidy(lm(FS ~ Days_Symptom_Onset_to_Visit_1, data = fs_pfs_df_with_manifest %>% dplyr::filter(parent == "4" & STIM == "Spike_2")))## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.0469 0.00600 7.83 1.68e-10
## 2 Days_Symptom_Onset_to_Visit_1 0.000129 0.000114 1.13 2.63e- 1# Note that the Healthies don't have values for Days_Symptom_Onset_to_Visit_1 and will get filtered out below
ggscatter(fs_pfs_df_with_manifest,
x = "Days_Symptom_Onset_to_Visit_1", y = "FS",
color = "STIM", palette = "jco",
add = "reg.line") +
facet_grid(parent ~ STIM) +
stat_cor()## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 46 rows containing non-finite values (stat_smooth).## Warning: Removed 46 rows containing non-finite values (stat_cor).## Warning: Removed 46 rows containing missing values (geom_point).
ggscatter(fs_pfs_df_with_manifest,
x = "Days_Symptom_Onset_to_Visit_1", y = "PFS",
color = "STIM", palette = "jco",
add = "reg.line") +
facet_grid(parent ~ STIM) +
stat_cor()## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 46 rows containing non-finite values (stat_smooth).## Warning: Removed 46 rows containing non-finite values (stat_cor).## Warning: Removed 46 rows containing missing values (geom_point).
For the CD4 FS vs Age plots, there is a positive association of Age with FS for Spike 1 (at least).
FS/PFS vs Sex
# Note that p-values on ggboxplot plots below are *not* adjusted, e.g.:
broom::tidy(wilcox.test(FS ~ Sex, data = fs_pfs_df_with_manifest %>% dplyr::filter(parent == "4" & STIM == "Spike_2")))## # A tibble: 1 x 4
## statistic p.value method alternative
## <dbl> <dbl> <chr> <chr>
## 1 369 0.239 Wilcoxon rank sum test two.sided# Note that future ggpubr update should make outlier.shape=NA occur by default when jitter is added, but include for now
ggboxplot(fs_pfs_df_with_manifest %>%
dplyr::filter(!is.na(Sex)),
x = "Sex", y = "FS",
color = "Sex", palette = "jco",
add = "jitter",
outlier.shape = NA) +
facet_grid(parent ~ STIM) +
stat_compare_means(aes(group = Sex))
ggboxplot(fs_pfs_df_with_manifest %>%
dplyr::filter(!is.na(Sex)),
x = "Sex", y = "PFS",
color = "Sex", palette = "jco",
add = "jitter",
outlier.shape = NA) +
facet_grid(parent ~ STIM) +
stat_compare_means(aes(group = Sex))
In many of the boxplots above, Males are shifted slightly higher on the FS/PFS axis compared to Females, but none of the comparisons are significant. Perhaps a t-test or un-stratified plot would be significant.
Hosp vs Non-Hosp
- Plot the COMPASS heatmaps again, this time including row annotations for Cohort
- Compare FS and PFS between Hosp and Non-Hosp
COMPASS Heatmaps
crList2 <- lapply(crList,
function(cr) {
cr$data$meta$Cohort <- factor(ifelse(cr$data$meta$Cohort %in%
c(NA, "Healthy control"), "Healthy", cr$data$meta$Cohort),
levels = rev(c("Healthy", "Non-hospitalized", "Hospitalized")))
cr
})
names(crList2) <- names(crList)
CohortColors <- c("Healthy" = "#dfc27d", "Non-hospitalized" = "#80cdc1", "Hospitalized" = "#018571")
row_ann_colors <- list(`Cohort` = CohortColors)
cytokine_annotation_colors <- rep("black", 30)plotCOMPASSResultStack(crList2[CD4RunNames],
row_annotation = c("Stim", "Cohort"),
variable = "Stim",
show_rownames = FALSE,
main = "CD4 COMPASS Runs",
fontsize = 15, fontsize_row = 15, fontsize_col = 13,
row_annotation_colors=row_ann_colors,
cytokine_annotation_colors=cytokine_annotation_colors,
order_by_max_functionality = T)## The 'threshold' filter has removed 6 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&IL17a Ax700&!IL2 PE&!IL4/5/13 APC&!TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&!TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&IL4/5/13 APC&!TNFa FITC
plotCOMPASSResultStack(crList2[NotCD4RunNames],
row_annotation = c("Stim", "Cohort"),
variable = "Stim",
show_rownames = FALSE,
main = "Not-CD4 COMPASS Runs",
fontsize = 15, fontsize_row = 15, fontsize_col = 13,
row_annotation_colors=row_ann_colors,
cytokine_annotation_colors=cytokine_annotation_colors,
order_by_max_functionality = T)## The 'threshold' filter has removed 2 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&IL17a Ax700&!IL2 PE&IL4/5/13 APC&!TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC
plotCOMPASSResultStack(crList2[c(CD4RunNames, NotCD4RunNames)],
row_annotation = c("Stim", "Cohort"),
variable = "Stim",
show_rownames = FALSE,
main = "CD4 and Not-CD4 COMPASS Runs",
fontsize = 15, fontsize_row = 15, fontsize_col = 13,
row_annotation_colors=row_ann_colors,
cytokine_annotation_colors=cytokine_annotation_colors,
order_by_max_functionality = T)## The 'threshold' filter has removed 3 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&IL4/5/13 APC&!TNFa FITC
I think a commonality to the subsets which are enriched in Hospitalized individuals is CD154.
Now look at each run individually
cytokine_order_for_annotation = c("IFNg V450", "IL2 PE", "TNFa FITC", "CD154 PE-Cy5", "CD107a PE-Cy7", "IL4/5/13 APC", "IL17a Ax700")
for(cr_name in c(CD4RunNames, NotCD4RunNames)) {
print(my_plot.COMPASSResult(crList2[[cr_name]],
y="Cohort",
#show_rownames = T,
fontsize=16, fontsize_row=16, fontsize_col=19,
row_annotation_colors=row_ann_colors,
cytokine_annotation_colors=cytokine_annotation_colors,
cytokine_height_multiplier=2.5, order_by_max_functionality = T,
main=cr_name,
cytokine_order_for_annotation=cytokine_order_for_annotation,
staircase_cytokine_annotation=T # Looks prettier and is consistent with dotplots
))
}## The 'threshold' filter has removed 7 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7## gTree[GRID.gTree.4657]## The 'threshold' filter has removed 4 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.4785]## The 'threshold' filter has removed 6 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.4922]## The 'threshold' filter has removed 10 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.5037]## The 'threshold' filter has removed 2 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.5164]## The 'threshold' filter has removed 3 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.5286]## The 'threshold' filter has removed 3 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7
## gTree[GRID.gTree.5417]## The 'threshold' filter has removed 6 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7
## gTree[GRID.gTree.5528]FS/PFS vs Cohort
# Note that p-values on ggboxplot plots below are *not* adjusted, e.g.:
broom::tidy(wilcox.test(FS ~ Cohort, data = fs_pfs_df_with_manifest %>% dplyr::filter(Cohort != "Healthy" & parent == "4" & STIM == "Spike_2")))## # A tibble: 1 x 4
## statistic p.value method alternative
## <dbl> <dbl> <chr> <chr>
## 1 162 0.000345 Wilcoxon rank sum test two.sided# Note that future ggpubr update should make outlier.shape=NA occur by default when jitter is added, but include for now
ggboxplot(fs_pfs_df_with_manifest %>%
dplyr::filter(Cohort != "Healthy"),
x = "Cohort", y = "FS",
color = "Cohort", palette = "jco",
add = "jitter",
outlier.shape = NA) +
facet_grid(parent ~ STIM) +
stat_compare_means(aes(group = Cohort)) +
scale_x_discrete(labels=c("Non-hospitalized" = "Conv\nNon-Hosp", "Hospitalized" = "Conv\nHosp"))
ggboxplot(fs_pfs_df_with_manifest %>%
dplyr::filter(Cohort != "Healthy"),
x = "Cohort", y = "PFS",
color = "Cohort", palette = "jco",
add = "jitter",
outlier.shape = NA) +
facet_grid(parent ~ STIM) +
stat_compare_means(aes(group = Cohort)) +
scale_x_discrete(labels=c("Non-hospitalized" = "Conv\nNon-Hosp", "Hospitalized" = "Conv\nHosp"))
Background-corrected magnitudes
crList.no_healthy <- lapply(names(crList2),
function(cr_name) {
cr <- crList2[[cr_name]]
print(sprintf("Accessing %s", cr_name))
cr$data$meta <- cr$data$meta %>%
dplyr::filter(!(`SAMPLE ID` %in% setdiff(cr$data$meta$`SAMPLE ID`, rownames(cr$fit$mean_gamma))))
stopifnot(all.equal(rownames(cr$fit$mean_gamma), cr$data$meta$`SAMPLE ID`))
stopifnot(all.equal(rownames(cr$data$n_s), rownames(cr$fit$mean_gamma)))
stopifnot(all.equal(rownames(cr$data$n_u), rownames(cr$fit$mean_gamma)))
stopifnot(all.equal(names(cr$data$counts_s), rownames(cr$fit$mean_gamma)))
stopifnot(all.equal(names(cr$data$counts_u), rownames(cr$fit$mean_gamma)))
not_healthy_idx <- which(cr$data$meta$Cohort != "Healthy")
cr$data$meta <- cr$data$meta[not_healthy_idx,] %>%
mutate(Cohort = factor(Cohort, levels = c("Non-hospitalized", "Hospitalized")))
cr$fit$mean_gamma <- cr$fit$mean_gamma[not_healthy_idx,]
cr$data$n_s <- cr$data$n_s[not_healthy_idx,]
cr$data$n_u <- cr$data$n_u[not_healthy_idx,]
cr$data$counts_s <- cr$data$counts_s[not_healthy_idx]
cr$data$counts_u <- cr$data$counts_u[not_healthy_idx]
cr
})## [1] "Accessing 4_VEMP"
## [1] "Accessing NOT4_VEMP"
## [1] "Accessing 4_Spike_1"
## [1] "Accessing NOT4_Spike_1"
## [1] "Accessing 4_Spike_2"
## [1] "Accessing NOT4_Spike_2"
## [1] "Accessing 4_NCAP"
## [1] "Accessing NOT4_NCAP"names(crList.no_healthy) <- names(crList2)
source(here::here("scripts/20200604_Helper_Functions.R"))
# Arial is used by dotplot function
# Arial font setup. Downloaded afms from https://github.com/microsoft/microsoft-r-open/tree/ec3fd89e5fb5794bd8149905c134ad801bb61800
Arial <- Type1Font(family = "Arial",
metrics = c(here::here("data/Arial_afm/ArialMT.afm"),
here::here("data/Arial_afm/ArialMT-Bold.afm"),
here::here("data/Arial_afm/ArialMT-Italic.afm"),
here::here("data/Arial_afm/ArialMT-BoldItalic.afm")))
pdfFonts(Arial = Arial)
# Note p-values are bonferroni adjusted in dotplots
bg_corr_dotplots <- pmap(.l = list(c(CD4RunNames, NotCD4RunNames),
rep(c("CD4+", "Not-CD4+"), each = 4),
list(c(-0.0015, 0.01), c(-0.002, 0.01), c(-0.003, 0.005), NULL,
c(-0.005, 0.007), c(-0.003, 0.008), c(-0.005, 0.013), NULL),
c(5, 5, 3, 5, 5, 5, 5, 5)),
.f = function(n, p, y, p_text_size) {
# "staircase effect" for categories df heatmap gets done automatically, unlike in my_plot.COMPASSResult
# Returned dotplot is cowplot
make_dotplot_for_COMPASS_run(crList.no_healthy[[n]], n, parentSubset = p,
return_output = T, current_ylim = y, include_0_line=T,
cytokine_order_for_annotation=cytokine_order_for_annotation,
p_text_size=p_text_size)
})## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"names(bg_corr_dotplots) <- c(CD4RunNames, NotCD4RunNames)
for(n in names(bg_corr_dotplots)) {
print(plot_grid(ggplot() +
theme_void() +
geom_text(aes(0,0,label=n), size=10) +
xlab(NULL),
bg_corr_dotplots[[n]]$Dotplot,
ncol=1, rel_heights = c(0.1, 1)))
}
sig_tests_tally <- lapply(names(bg_corr_dotplots), function(n) {length(which(bg_corr_dotplots[[n]]$Test_Results$p.adj < 0.05))})
names(sig_tests_tally) <- names(bg_corr_dotplots)
sig_tests_tally # make sure all significant comparisons were plotted and accounted for## $`4_Spike_1`
## [1] 8
##
## $`4_Spike_2`
## [1] 4
##
## $`4_NCAP`
## [1] 5
##
## $`4_VEMP`
## [1] 0
##
## $NOT4_Spike_1
## [1] 0
##
## $NOT4_Spike_2
## [1] 0
##
## $NOT4_NCAP
## [1] 0
##
## $NOT4_VEMP
## [1] 0Focus on IL2
Look at each run individually again, this time with IL2 subsets on RHS (since there was evidence of elevated DMSO IL2 signal in hosp)
for(cr_name in c(CD4RunNames, NotCD4RunNames)) {
print(my_plot.COMPASSResult(crList2[[cr_name]],
y="Cohort",
#show_rownames = T,
fontsize=16, fontsize_row=16, fontsize_col=19,
row_annotation_colors=row_ann_colors,
cytokine_annotation_colors=cytokine_annotation_colors,
cytokine_height_multiplier=2.5, order_by_max_functionality = T,
cytokine_order_for_annotation = c("IL4/5/13 APC", "IFNg V450", "TNFa FITC", "IL17a Ax700",
"CD154 PE-Cy5", "CD107a PE-Cy7", "IL2 PE"),
dichotomize_by_cytokine = "IL2 PE",
staircase_cytokine_annotation = TRUE,
main=cr_name))
}## The 'threshold' filter has removed 7 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7## gTree[GRID.gTree.8731]## The 'threshold' filter has removed 4 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.8859]## The 'threshold' filter has removed 6 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.8996]## The 'threshold' filter has removed 10 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.9111]## The 'threshold' filter has removed 2 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.9238]## The 'threshold' filter has removed 3 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.9360]## The 'threshold' filter has removed 3 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7
## gTree[GRID.gTree.9491]## The 'threshold' filter has removed 6 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7
## gTree[GRID.gTree.9602]Previously, we observed higher DMSO IL2 signal in Conv Hosp compared to Conv Non-hosp, but so far there’s no obvious sign that this difference carries over universally to the background-corrected IL2-containing subsets. Quantitative magnitude boxplots should provide more detail.